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ABSTRACT 

Experimentally obtained dynamics of time-optimal, horizontal head 
rotations have previously been simulated by a sixth order, non- 
linear model driven by rectangular control signals. EMG 
recordings have aspects which differ in detail from the 
theoretical rectangular pulsed control signal. We have obtained 
control signals for time-optimal as well as sub-optimal 
horizontal head rotations by means of a newly developed inverse 
modelling proceedure. With experimentally measured dynamical 
data serving as the input, this procedure inverts the model to 
produce the neurological control signals driving muscles and 
plant. The relationships between these controller signals, and 
EMG records should contribute to our understanding Of the 
neurological control of movements. 
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INTRODUCTION 


Head movements are similar to arm movements about the elbow in 
dynamics and time scales (Lehman 1983) and are of interest 
because of their interactions with eye movements, posture, and 
perception. Zangemeister , et. al., 198la-e, have studied head 
movements and their involvement in shifts of gaze, ’the eye's 
position in space. They have also quantified the dynamics of 
time-optimal horizontal head rotations in terms of the peaks of 
the dynamical variables position, velocity, and acceleration and 
plotted them in the Main Sequence diagram to show the 
relationship between dynamics and movement magnitude 


Interest in the control mechanisms involved in head movements has 
lead to the study of the electromyographic activity of neck 
muscles involved in these movements (Zangemeister, Stark, 
Meienberg, & Waite & Stark, Hannaford et al. 1983) and -to efforts 
to model the system. 


* On Leave: Department of Electronics Engineering, Kwang Woon 
University, Seoul, Korea 
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Figure 1 

6th Order Non-Linear Model of the head rotation system. 
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Zangeraeister Lehman and Stark (8 lab) simulated the muscles and 
plant of the head movement system with a 6th order non-linear 
model incorporating Hill’s force-velocity relationship, two 
antagonist muscles, and an overdamped second order plant (Figure 
1). Their model matched experimentally measured Main Sequence 
dynamic peaks when driven by heuristically derived control 
signals. Versions of this model have had a fruitful history of 
application to many different physiological systems. Stark 
(1961) proposed and Atwood et. al. (1961) simulated a two-muscle 
model for understanding neurological control mechanisms. Cook & 
Stark (1968), and Clark & Stark (1975) used more detailed 
versions with appropriate parameter values to model saccadic and 
other eye movements, and it has been applied to the eyelid in 
modeling the dynamics of the blink (Kim, et. al. 1984b). 

In some cases, it is possible to invert a numerical model and 
obtain controller signals as a function of dynamics. Cook (1965) 
linearized the model and obtained a closed form solution to the 
inverse problem. Recently, Kim, et. al.(84a) has developed an 
iterative method which has been applied to the non-linear model 
of the eye-movement system. An adaptation of this technique was 
used in the present work. 

METHODS 

Experimental 

Horizontal head position was measured by a precision 
potentiometer attached to a bicycle helmet frame worn by the 
subject. To allow for vertical head movement and subjects of 
differing heights, the potentiometer was coupled to the head 
through a compliant fitting which unfortunately resulted in a 
delay of about 50 ms. between actual and recorded head position. 


Electromyographic activity was measured differentially with two 
S&W number 737 self-adhesive Ag-AgCl disposable electrodes placed 
approximately 5cm apart on the skin along the major axes of the 
left and right splenius capitus muscles. EMG and position data 
was digitized at 1 000 Hz by an LSI 11-23 computer with 12 bit 
analog to digital converters (Hannaford et.al. 1983). 


The subject's head movements were made in response to light 
emitting diode (LED) targets alternately flashing at points on a 
curved screen 1 meter from the subject's head. The subject was 
aware of the exact position of his head by a small spot of light 
projected from his helmet onto the screen. When the target 
illumination alternated between the two positions, (at intervals 
of 4 seconds) the subject performed 20, 40, and 60 degree 
horizontal head movement’s. 


The subjects were instructed to move their heads "as fast as 
possible" to produce an intent to respond to the target in a 
time-optimal manner. This experiment resulted in stereotyped 
movements which could be ensemble averaged along with their 
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rectified EMGs . 


In a separate experiment, the targets were set 40 degrees apart 
and the subject was given different instructions for movements in 
the two different directions. For movements to the right, the 
subject was asked to make time-optimal movements as in the first 
experiments. For the leftward movements, he was instructed to 
move "however you want." Of the many and varied leftward 
movements that resulted from this paradigm, two leftward 
movements occurring directly after each other (with one time- 
optimal rightward movement between them) were selected for 
analysis . 

Modeling 

For simulation of the horizontal head rotation system, we used 
the sixth-order non-linear model developed by Zangemeister , 
Lehman, and Stark ( 1 98lab). This model consists of two 
identical, antagonistic muscle elements driving a second order 
plant (Figure 1). The muscle elements have a force generator 
driven through a first order low pass filter representing the 
calcium activation process. The control signals, nl and nr, 
range from zero to one to represent the possible range of 
excitation from none to full excitation. To help the reader's 
intuition, we have plotted this signal in equivalent kilograms to 
suggest the steady state force that would result from constant 
excitation at a given level. In parallel with the force 
generator is a non-linear viscous element representing the Hill's 
force-velocity relationship (Hill, 1938). Force is transmitted 
to the load through a series elastic element representing the 
properties of muscle tendon and attached cross-bridges. 

This system is modeled and simulated by a set of 6 state 
equations and two ancillary equations (Table 1). The values used 
for the parameters (Table 2) are based on previous work 
(Zangemeister, Lehman & Stark, 1981a) and recent, improved 
estimates (J.M. Winters, private communication). 
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State Equations 


x = v 

vl = (htl - ks(xl-x)) / bl 
vr = ( -htr- ks(xr-x)) / br 

a = (-kp x - bp v + ks (xl - x) + ks(xr - x) ) / j 
dhtl = (nl-htl)/ta 
dhtr = (nr-htr)/ta 

Ancillary Equations 

bl = (1.25 htl) / (bh + vl) 

1.25 htl / 900. 


br = (1.25 htr) / (bh + vr) 

1.25 htr / 900. 


TABLE 1 

Equations for the sixth order non-linear model 
of horizontal head rotation. 


vl > 0 
vl < 0 


vr < 0 
vl > 0 


Name 

Symbol 

Value 

Units 

Parallel Viscosity 

bp 

2.0 

gr-f deg _1 sec 

Activation Time Const. 

ta 

50.0 

milli seconds 

Hill's Constant; b 

bh 

350.0 

deg sec -1 
gr-f deg - }sec 
gr-f deg“' 
gr-f deg” 1 

Rotational Inertia 

j 

0.18 

Series Elasticity 

ks 

350.0 

Parallel Elasticity 

kp 

2.0 


TABLE 2 



Model 

Parameter 

Values. 



Inversion bv Iteration of Forward Model 

The simulations typically undertaken with such a model involve 
applying various control signals at the model's inputs (the 
neurological force commands to the muscles, nl and nr, are 
presumably a product of firing rate and recruitment) and 
observing the model's output; position, velocity, and 
acceleration time functions. In this study, we reversed the 
process, obtaining the control signal as a function of 
experimentally recorded movement dynamics. 

To do this, we used an iterative method which, at each time 
sample, involves finding the control signal values which result 
in model output that exactly matches the experimentally recorded 
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Figure 2 

Block diagram of iterative model inversion process. 
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dynamics at that time (Figure 2). When these 

values are obtained, the state variables are updated and the 
IJoJessis repeated for the next time sample. Several issues 
arise which must be resolved before this computation can be 

performed . 

Because the system has only one measurable output, head position, 
only one independent variable can be obtained by inversion of 
model In order to use an iterative method to minimize the 
difference between model output and experimental data, the two 
model inputs, agonist and antagonist control signals must be 
constrained to be a function of this single variable. In these 
simulations, the independent variable was net force (fnet) and 
the constraint used was: 


-.agonist 

Antagonist 

Agonist 

Antagonist 


fmin 

fmin 

fmin 

fmin 


fnet 


fnet 


fnet > 0 


fnet <= 0 


suggested by the absence of co-contraction shown m the EMG 
recordings we analyzed. 

Iteration Methods 

with net force driving the model through the constraints and 
generating' Agonist? ^antagonist force commands, the problem 
becomes to find the value of net force for which 

E = 0 


where 

E = V h (t) - V m (t,fnet) 

and V h (t) is head velocity as a function of time over (0 < t < 
tmax), h and V m (t,fnet) is the model output as a a function of fnet 
and its previous values. Solving this equation for each value of 

t from 0 to tmax yields net force as a function of time. y 

applying the constraint, we then have agonist and antagonist 
force control signal as a function of time. 

The method initially used to solve this equation at each time 
sample was an adaptation of the Newton-Raphson method m which 
fnet was iterated until the value of E was less than * 

epsilon. After E was calculated for two values of fnet, the new 

estimate for fnet is 
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where j is the iteration number. 


Although this method was effective, it occasionally failed to 
converge when V m (t,fnet) was sufficiently non-linear as a 
function of fnet. It was subsequently found that a binary search 
method would guarantee convergence of the algorithm. 

In this second method, an initial range of values is selected 
between which it is assumed must lie the correct value of fnet. 
This range can easily be determined by taking the maximum 
expected force value and allowing fnet to vary between that value 
both above and below zero. The first estimate in this procedure 
is zero. Then each subsequent estimate is improved by an 
increment equal to the maximum value divided by a larger and 
larger power of two. If the value of E resulting from this new 
estimate of fnet is negative, the next increment is subtracted 
from fnet. If it is positive, it is added. Convergence relies 
on the assumption that E crosses zero at least once at some value 
of fnet between the initial guesses, an assumption which can 
always be made true by widening the initial range at a slight 
expense in convergence time. 

Numerical precision 

Compared to the eye movement system (Kim et al, 1984a), the head 
movement system has very long time constants; a step change in 
controller signal results in a very small instantaneous change in 
head velocity. Also, small amounts of noise, including 
quantization noise, in the velocity signal will require large 
changes in the control signal in order for the model output to 
match this noise. Thus attention must be paid to numerical 
precision and filtering if this calculation is to be successful. 
In our computations, we used double precision arithmetic for all 
calculations of state variables, ancillary equations, and system 
error, E. Furthermore, we filtered the input data to produce a 
double precision result with a sufficiently small amount of noise 
and quantization error. Filtering did not result in appreciable 
changes in movement dynamics. 

Filtering and Effects of Bandwidth 

Filtering of head movement trajectories was necessary to reduce 
undesired measurement and quantization noise. For this smoothing 
operation, a Hamming window, zero-phase-delay, low pass filter 
was used (Rabiner & Gold, 1 975). The frequency response of the 
ideal low pass filter is 

H(jw) = 1 for -wc < w < wc 
0 elsewhere 

where wc is the desired cutoff frequency. The filter we used is 
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Figure 3 

Smoothing of movement dynamics necessary for model inversion was 
performed with a 100 point Hamming window low pass filter. 
Impulse response of this filter (top) is symmetrical about zero 
to eliminate phase delay. Frequency response (bottom) shows 
sharp cutoff at 20 Hz. 






one which comes close to this ideal response. The impulse 
response of the ideal low pass filter is 


n wc 

sin — . • .i — - — 
ws 

h(n) = — - 

pi n 

where ws is the sampling frequency. The ideal impulse response 
extends to infinite values of n. The Hamming window is employed 
as a finite weighting sequence on the infinite ideal impulse 
response to produce smooth truncation. The weighting function of 
the Hamming window is 

w(n) = 0.54 + 0.46 cos(2pi n/ N), -N n <= N 

where N is the number of data points for the truncation window. 
The modified impulse response weighted by the Hamming window is 

hw( n) = h( n) w( n) . 

The output sequence y(n) of the Hamming window low pass filter is 
given by the convolution of the input sequence with the modified 
impulse response hw(n). Note that hw(n) is symmetrical with 
respect to hw(0). The filtered output sequence y(n) can thus be 
described by a finite difference equation as; 

y(n) = hw(0) x(n) + hw(1) (x(n-1 )+x(n+1 ) ) 

+ hw(2) (x(n-2)+x(n+2) ) + ... 

For smoothing the position and velocity trajectories, we used a 
Hamming window low pass filter of 100 data point with a cutoff 
frequency of 20 Hz at a sampling rate of 1000 Hz (Figure 3)- 


RESULTS 


Time Optimal Movements 

We prepared 3 ensemble averages of fast movements at amplitudes 
of 20, 40 and 60 degrees (n = 5). Velocity and acceleration 
traces show amplitude-dependent peak values characteristic of 
time-optimal movements (Figure 4). Full-wave-rectified EMG 
activity from agonist and antagonist muscles was also averaged 
(Figure 5). The EMGs exhibit the tri-phasic burst pattern found 
in fast movements about several different joints (Wachtolder, 
Angel, Ghez & Martin, Litvintsev & Seropyan, Wadman, van der Gon, 
& Derksen, Hannaford et. al., 83, Cheron & Godaux). It is 
difficult to quantify signals of this type in terms of height and 
width. However, as a function of movement magnitude, they seem 
to vary more in width than in amplitude although PA, the first 
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Acceleration (deg/sec**2) Velocity (deg/sec) Position (deg) 



Figure 4 

Three ensemble averages of time-optimal horizontal head movements 
of 20, 40, and 60 deg. (n=5). Also shown are velocity (middle) 
and acceleration (bottom) computed from the filtered data. 
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Mm« r ^ i: t e - r ® ignais obtained by inverse modelling of the three 
°P fclmal movements of figures 4 and 5 . Three pulses of 
excitation are seen which correspond to the three EMG pulses 
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agqnist EMG burst reaches a peak value about twice that of PB and 
PC. 

The result of the inverse modelling process is a pair of control 
signals describing the excitation levels of the antagonistic 
muscle pair for each of the three movement magnitudes. This 
signal (Figure 6) also shows the tri-phasic pattern, having an 
initial agonist burst followed by a burst of antagonist activity 
and finally a second antagonist burst. PA, the first agonist 
pulse, increases from 75 to 100 ms in width and from 6 to 17 
equivalent kilograms (kge) in amplitude as movement magnitude 
increases from 20 to 60 degrees. PB ranges from 70 to 100 ms in 
width and has a relatively constant amplitude of about 6kge 
except for a second peak of about 8 kge in the 60 degree case. 
PC appears to be fairly constant at about 40 ms in width and 
about 2 kge in amplitude. The skew evident in PB and PC is due to 
the concatenation of the three pulses; as longer pulses are 
concatenated, later pulses are delayed appropriately. 

Comparison of the EMG records with controller signals resulting 
from the inversion shows a delay of about 30 ms. between EMG and 
control signal resulting from delay in the head position 
measurement apparatus not accounted for in the model. The width 
of PA, the first agonist EMG burst matches well with the first 
agonist control signal pulse for all three magnitudes. The 
antagonist excitation pulse shows the same increase in onset 
times (of roughly 10 ms./20deg) with movement magnitude as does 

PB, the antagonist EMG burst. But each antagonist control signal 
pulse is longer than the corresponding EMG pulse. While the EMG 
pulse, PA is of roughly constant amplitude, the first control 
signal pulse amplitude varies over a three to one range with 
movement magnitude. 

PC, the second agonist EMG burst, is of roughly constant 
amplitude but its width varies strongly with movement magnitude 
from about 60 to 75 ms. Width of the third control signal pulse 
is difficult to ascertain because of its approximately 
exponential decay, but unlike the PC of the EMG, its amplitude is 
quite small relative to the first pulse at all three magnitudes. 
A hypothetical linear relationship between EMG magnitude and 
control signal magnitude would suggest that the third control 
signal pulse have an amplitude roughly half of that of PA but in 
fact, it is much less. This may suggest that the model is 
too viscous during the later phase of the movement, requiring 
less active damping than does the real system. 

A Pair of Slow Movements 

Figure seven depicts two double movements. These single records 
(not averages) were taken of the subject on successive leftward 
rotations under instructions to move "however you want to". Both 
records begin with a movement of approximately 15 degrees to the 
left which is of insufficient amplitude to reach the target. The 
peak velocity of this movement was 210 deg/sec for the solid and 
220 deg/sec for the dashed record. A second deflection of the 
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Antagonist EMG CuV) 



Time (ms) 


Figure 7 

Records of two double movements and their EMG's. First movement 
is nearly identical in both records while a second, corrective, 
movement appears in the later record (solid trace) about 240 ms. 
after the first. Note similarity in duration, amplitude and 
number of pulses of EMG bursts causing first movement, and 
divergence of EMG’s as dynamics diverge. 
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position trace appears about 100 ms. later. About 220 ms after 
the first movement, the subject made a second, corrective 
movement to cover the remaining distance. The solid record makes 
a discrete corrective movement of about 15 deg. with a velocity 


peak of 135 deg/sec. In the dashed trace, there is instead a 
drift at an approximately constant velocity of about 20 deg/sec. 
The first movement is nearly the same in both records while the 
second movement is of larger amplitude and velocity in the second 
of the two records (solid trace). 


An initial burst of agonist EMG is seen in both records. These 
bursts are similar in amplitude, duration, and number of spikes. 
The antagonist channel shows a small amount of EMG activity due 
to either cross-talk or a small amount of co-contraction. No 
antagonist EMG activity appears after the initial agonist burst 
(PA). The second EMG burst is much more prominent in the second 
record (solid) and the corresponding second movement is greater. 


Control signals were calculated using the inverse model on the 
velocity trajectories of the two records (Figure 8). The 
calculated control signals consist of a series of rounded, 
roughly triangular pulses, the first and largest ones resulting 
from the initial movement in both records. These pulses have a 
peak force value of about 5 kge. with the dashed pulse slightly 
greater coresponding to the slightly greater peak velocity. Both 
pulses are about 60 ms in duration. The slight difference in 
peak force corresponds to the slightly faster time course of the 
first movement in the first record (dashed lines). The control 
signal obtained by the inversion contains 6 subsequent, smaller, 
pulses of activity in the agonist and 6 in the antagonist. 


The second agonist pulse (about 1 kge peak force, and about 50 ms 
duration) corresponds to the slight increase in velocity seen 
about 100 ms. before the start the second movement. 220 ms after 
the initial agonist burst, the second movement is initiated by 
another burst of agonist activity. Here, as the dynamics diverge 
between the two records, a larger agonist burst (3 kge vs less 
than 1 kge peak force and 60 ms vs 50 ms pulse width) appears in 
the later (solid) record corresponding to the larger movement. 


Small pulses of antagonist force follow immediately after each 
agonist pulse and immediately precede the next agonist pulse. 
These pulses are not present in the EMG records and their 
presence may be due to the fact that the filtering of the 
dynamics results in slight dynamical changes requiring a smoothly 
changing or "ringing” control signal. Another possibility is 
that since no co-contraction is possible because of the 
control signal constraint, active damping may be required of the 
control signal to make up for lack of an active viscosity due to 
co-contraction but no co-contraction is evident in the EMG. 
Finally, it may be that other neck muscles such as the left and 
right sterno cleido mastoid may be activated at these times. 
Recording of additional EMG channels may clarify this possiblity. 
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Figure 8 


Control signals obtained by inverse modelling of two slow 
movements shown in figure 7. Timing shift between EMG's (fig. 7 ) 
and control signal pulses (this fig.) is due to time delay in 
experimental apparatus. Extra pulses of controller signal 

nnnVlnt ^ m * y b -, e d u e ^ 0 la ° k of c o - c o n t r a c t i o n imposed by 
controller signal constraint, or activity of other neck y 

not recorded in fig. 7 . 617 


muscles, 





DISCUSSION 


Inversion of the head model is an interesting problem in 
numerical analysis. Because of the long time lags involved, the 
problem is near to being ill conditioned (Rice, 1 984). Use of 
FORTRAN'S double precision arithmetic was required both for the 
computation of state variables and for the results of the data 
filtering against which model output was matched. 

Algorithms for solving for a zero of the output error function E 
must be sufficiently able to deal with the non-linear behavior of 
the model to guarantee convergence in a reasonable amount of 
time. A binary search algorithm was found to always converge but 
to take more time, in many cases, than a successive approximation 
method based on the Ne wton-Raphson algorithm. A suitable 
improvement would be a successive approximation method generating 
a rough estimated range, followed by a binary search to guarantee 
convergence and full double precision accuracy. 

Another area for further study is that of the controller signal 
constraints by which net force is converted to agonist and 
antagonist control signals. The constraint used in this study 
was suggested by the fact that appreciable co-contraction was not 
evident in the EMG records from time-optimal movements. Cook 
(1965) and Kim, et.al. (1984a) in inverting the eye movement 
system have used another possible constraint, which specifies a 
given amount of co-contraction in terms of a ratio of antagonist 
to agonist excitation. For example, if the co-activation level 
is set to 20%, and the excitation level required is 1, then the 
agonist would be set to 1 and the antagonist to 0.2. It will be 
interesting to recompute the above results with this type of 
constraint and assess the effect of co-contraction level on 
antagonist activity in the slow movements. 

The inverse modeling process is an aid to understanding the 
control of skeletal muscle in movements and helps create a 
conceptual link between EMG and movement dynamics. Completion of 
this link will yield the ability to predict movement dynamics 
from a knowledge of the plant and of the EMG signal. The three 
unresolved steps in this link are the calibration of EMG to 
excitation; the improvement of experimental apparatus to reduce 
delays, non-linearities, and frequency dependent effects; and the 
further elaboration of the subtle non-linearities in the model. 

The above process of "Dynamic Calibration" would be a step 
toward an ideal EMG signal processor in the sense that the 
response of the Hamming window low pass filter approaches that of 
the ideal low pass filter. So far, the only signal processor to 
adequately interpret the EMG is still the living muscle. 
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